Grid-based simulation program for gravitational wave interferometers 
with realistically imperfect optics 
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We describe an optical simulation program that models a complete, coupled-cavity interferometer 
like those used by the Laser Interferometer Gravitational- Wave Observatory (LIGO) Project. A wide 
variety of interferometer deformations can be modeled, including general surface roughness and sub- 
strate inhomogeneities, with no a priori symmetry assumptions about the nature of interferometer 
imperfections. Several important interferometer parameters are optimized automatically to achieve 
the best possible sensitivity for each new set of perturbed mirrors. The simulation output data set 
includes the circulating powers and electric fields at various points in the interferometer, both for 
the main carrier beam and for its signal-sideband auxiliary beams, allowing an explicit calculation of 
the shot-noise- limited gravitational- wave sensitivity of the interferometric detector to be performed. 
Here we present an overview of the physics simulated by the program, and demonstrate its use 
with a series of runs showing the degradation of LIGO performance caused by realistically-deformed 
mirror profiles. We then estimate the effect of this performance degradation upon the detectabil- 
ity of astrophysical sources of gravitational waves. We conclude by describing applications of the 
simulation program to LIGO research and development efforts. 
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I. INTRODUCTION 

The LIGO Project is part of the current initiative 
to detect gravitational radiation via its perturbing effects 
on resonant laser interferometers. This initiative involves 
several large collaborations around the world, including 
VIRGO 0, GEO 0, TAMA 0, and ACIGA Q. 

A paramount issue in all such projects is the maximiza- 
tion of interferometer sensitivity to detect the extremely 
weak signals that are expected from even the most power- 
ful astrophysical sources To that end, coupled-cavity 
systems with multiple resonant stages will be used to 
maximize the shot-noise-limited signal-to-noise ratio of 
the gravitational wave (GW) signal readout. For these 
complex interferometric detectors, intensive modeling is 
necessary to estimate their performance in the presence 
of general optical imperfections. 

A hierarchy of approaches has been used to esti- 
mate detector performance, each method negotiating the 
trade-off between accuracy and computational complex- 
ity. Analytical methods may suffice for the consideration 
of optical defects that can be treated as pure losses, or 
for some cases involving geometric or randomized .§j 
mirror deformations. A matrix model that evaluates the 
coupling between the first few lowest-order TEM laser 
modes has been shown to be useful for the study of mir- 
ror tilts and beam displacements Q; and another ma- 
trix model using discrete Hankel transforms exists for 
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problems with axial symmetry |lO|| . Such models, which 
consider the exchange of power between a limited set of 
pre-specified modes, allow one to obtain fast results that 
can be used for predicting certain important interferom- 
eter behaviors (e.g., time-dependent detector responses 
in in El), at the expense of some sophistication in 
modeling the detailed interferometer steady-state power 
buildup. For the consideration of highly general optical 
imperfections, however, the most comprehensive method 
is the complete modeling of the transverse structure of 
the laser field wavefronts; this method simulates the elec- 
tric fields on large grids, and uses (Fourier transform- 
based) numerical computations for the propagations of 
these laser beams through long cavities [l4L Il5l |. This 
technique is useful in the case of mirrors with complex de- 
formations, and for mirrors with significant losses (due to 
absorption, scattering, and/or diffractive loss from finite- 
sized apertures), that violate the assumption of unitarity 

Pi 

for mirror operators in the matrix models, thus intro- 
ducing complications into that approach ^3- 

Grid-based simulations of intra-cavity laser fields have 
a long history (e.g.. |l6| ) , and have been applied previ- 
ously (e.g., [l7lllMll9| ) to the study of the interferometric 
GW detectors now being implemented. But simplifica- 
tions are generally imposed, such as restricting the opti- 
cal deformations that are studied (e.g., considering only 
geometric imperfections, like tilt and curvature mismatch 
,17, 18)], and/or by modeling a relatively simple cavity 
system [T9j. In this paper, we describe a simulation pro- 
gram (originally based upon the work of Vinet et al. 'l5j) 
that has been extended to efficiently model the complex 
^^ Ids that build up between realistically imperfect optical 
components, while resonating within complete, coupled- 



cavity interferometers like those used for the LIGO detec- 
tors. Versions of our program have been used for a vari- 
ety of appHcations by the gravitational wave community, 
including numerous design and performance estimation 
tasks conducted by the LIGO group itself (see Section 
IVIfl . as well as for collaborative investigations between 
LIGO scientists and other groups — such as ACIGA- 
LIGO efforts to explore alternative interferometer length 
control schemes and TAMA-LIGO efforts to esti- 
mate the effects of mirror imperfections and thermal 
lensing upon the performances of their future, large- 
scale interferometers. 

A wide variety of interferometer imperfections can be 
modeled with our program, including mirror tilts and 
shifts, beam mismatch and/or misalignment, diffractive 
loss from finite mirror apertures, mirror surface figure 
and substrate inhomogeneity profiles — in particular, us- 
ing deformation phase maps that are adapted from mea- 
surements of real mirror surfaces and substrates — as 
well as fluctuations of reflection and transmission inten- 
sity across the mirror profiles. In addition to the main 
(carrier frequency) laser field, auxiliary fields (i.e., radio 
frequency sidebands) for LIGO's heterodyne signal de- 
tection scheme are also modeled, allowing us to give ab- 
solute numbers for the shot-noise-limited GW-sensitivity 
of a single LIGO interferometer. Furthermore, a number 
of active optimization procedures are performed contin- 
uously during code execution, guaranteeing that several 
key parameters in the LIGO configuration will be brought 
to their optimum values upon program completion; the 
GW-sensitivity is thus maximized for each specific "real- 
istic" interferometer, given the particular imperfections 
being simulated in that run. Considering the time needed 
for program execution, we note that though our prin- 
cipal usage of the program has been on various super- 
computing platforms, a full run of the simulation code 
(with all options included) can be performed in a non- 
prohibitive amount of time on a modest Sun SPARCsta- 
tion, even for significantly (and non-symmetrically) de- 
formed mirrors. This has been achieved through a com- 
bination of fast iteration techniques, efficient parameter 
optimization routines, and procedures designed to care- 
fully choose the initial guesses for laser fields that are 
computed via relaxation. Lastly, we note that versions of 
the code are available for both the first-generation LIGO 
and advanced-LIGO (Dual Recycling "2^) configurations 
E3i l25j , though we will focus primarily upon the former 
in this paper. 

The discussion is organized as follows: in Section ^1 
we give a description of the physical system that is mod- 
eled here (i.e., a first-generation LIGO interferometer), 
and the computed electric fields that are necessary for 
the calculation of its shot-noise-limited GW-sensitivity 
function. In Section IIIII we provide an overview of the 
technical details of the program's operations, including 
the modeling of the optics, the iterative method used 
for computing the interferometer's electric fields, and the 
various optimization procedures that are performed to 
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FIG. 1: Schematic diagram of the core optical configuration of 
a first-generation LIGO interferometer. (Not drawn to scale.) 



maximize its sensitivity. In Section IIVI we present an 
array of runs with representative sets of mirror substrate 
and surface deformation maps adapted from real mirror 
measurements, in order to demonstrate how realistically- 
deformed mirrors will reduce the circulating power stored 
in a LIGO interferometer (as well as altering its modal 
structure), thus degrading the interferometer sensitivity 
and interfering with its control systems. In Section^ we 
discuss the impact of optical imperfections upon LIGO 
science capabilities, by estimating how deformed-mirror 
effects will reduce LIGO's ability to detect astrophysical 
sources of gravitational waves, such as non-axisymmetric 
pulsars and coalescing black hole binaries 6]. In Sec- 
tion ^1 we conclude with a discussion of various LIGO 
research and development initiatives to which this grid- 
based simulation work has contributed. 



II. THE PHYSICAL SYSTEM THAT IS 
MODELED 

Figure ^ is a schematic diagram of the core optical 
configuration of a full-sized, first-generation LIGO inter- 
ferometer 26] . A summary of typical interferometer (and 
computational) parameters for this configuration is pre- 
sented in Table m These are the primary program input 
values that we will use for the sample runs to be pre- 
sented in Section HVI Not shown (or modeled by us) are 
the mode cleaning, frequency stabilizing, and matching 
optics which prepare the laser light for the interferome- 
ter; also not modeled are the pickoffs, phase modulators, 
or the full control system apparatus that will be used in 
a real interferometer to read out its complete operational 
state prj . 

The system depicted in Fig. essentially a Michel- 
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TABLE I: Typical parameter values for a first-generation 
LIGO interferometer, including physical specifications and 
computational parameters required for the simulation pro- 
gram. The labeling of the optical elements are as depicted in 
Fig-Q Some parameters are optimized during code execution, 
and are given here only as approximate ranges of values. 



Quantity 


Value(s) 


Laser Wavelength 


1.064 ^m (Nd:YAG light) 


Modulation Frequency 


u^od ~ 24.0 MHz 




Li = 5.0 m 




L2 = 4.19 m -1- Lasymm 


Cavity Lengths 


L3 = 4.19 m - Lasymm 



iasymm ~ 9 — 25 Cm 

L4 = L5 = 4.0 km 



-^curv, 1 

= 10.0 km 

Mirror Curvature Radii -Rcurv,2 = i?curv,3 = 14.6 km 

Rcuvv,4, = ^curv,5 — 7.4 km 

i?i ~ .9861 - .9390 
Mirror Intensity Reflectivities R2 ~ Ra = .97 
(Refl.-Side) R4 = R5 ^ .99994 

Rbs = .49992 

Mirror Intensity Reflectivities Ri ~ .9861 - .9390 
(AR-Side) R2 = R3 ^ .968817 

Rbs = -49971 

Mirror Intensity Transmissions Ti ~ .01385 - .06095 
(Both Sides) T2 = T3 = .02995 

(Pure Loss =1- R-T) Tbs = .50003 

Beam Waist Diameter 7.0 cm 

24 cm (Circular Profiles), 
Mirror Aperture Diameters 24.4 x 17.2 cm (Beamsplitter 

at 45° w.r.t. Beam Axis) 
Mirror Thicknesses Beamsplitter = 4 cm 

(Perpendicular to Surface) All Others = 10 cm 
Substrate Refraction Index n = 1.44963 
Calculational Window Size 70 cm x 70 cm (Square) 
Gridding of Calc. Window 256 x 256 pixels 



son interferometer, converts the differential arm length 
changes caused by a gravitational wave (GW) into an os- 
cillating output field amplitude at the exit port of the 
beamsplitter, where a carrier field dark-fringe would oth- 
erwise (ideally) have been maintained. The partially- 
transmitting input mirrors (r2,73) and highly-reflective 
end mirrors Rq) form Fabry-Perot (FP) arm cav- 
ities in the "inline" and "offline" arms, which amplify 
this effect; and the power recycling mirror (Ri) provides 
a broadband ampliflcation for the stored energy in the 
interferometer that is available for signal detection [28l| . 
Coupling occurs between the fields from each of the dif- 
ferent cavities in this optical configuration, and is respon- 
sible for the complex response behaviors of the system. 

Radio- frequency modulation sidebands ("RF- 
sidebands" ) impressed upon the carrier-frequency 
input light serve as a "local oscillator" for a heterodyne 
detection scheme |23| • Substantial RF-sideband power is 
made to emerge from the beamsplitter exit/signal port, 
and it is added to the quantity of GW-induced carrier 
light that is escaping there, to ultimately produce an 



output signal linear in /i, the dimensionless amplitude 
of the GW. The carrier power is maximized by giving 
it a double-resonance: it is first resonant in (either of) 
the FP arm cavities, and again resonant in the "power 
recycling cavity" (PRC) formed by mirrors R2, and 
_R3. Alternatively, the RF-sidebands (in this standard 
LIGO detection scheme) are only resonant in the PRC, 
to take advantage of the broadband amplification there; 
but they are far-off-resonance in the long FP-arms, 
so that they may serve as a stable reference which is 
unaffected by gravitational waves. 

Assuming that the carrier is held to its double reso- 
nance (giving it a phase shift of tt in resonant reflection 
from the FP-arms |2a|), the RF-sidebands have their own 
resonance requirements, defined by the two conditions: 
(2 fcmod iarm) « A^odd ' 71", and (2 k^od Lpjic) ~ A^odd • TT. 
Referring to Fig. Larm represents L4 or L5, Lp^c = 
Li + {L2 + Ls)/2, and the RF-modulation frequency is de- 
fined according to k^od = STrfmod/c = 27r|i/carr - i^sb\/c. 
The first (anti-resonance) condition is not enforced as a 
precise equality, in order to avoid the unwanted resonance 
of non-negligible second-order modulation sidebands (at 
2 i^mod) in the FP arm cavities. But the second (PRC- 
resonance) condition must be achieved to sub- wavelength 
tolerances, and requires careful fine-tuning corrections 
due to the effects of realistically-deformed optics (see 
Sec. liii(J2ll . 

The heterodyne GW-signal is obtained by interfer- 
ing the GW-induced carrier beam at the beamsplitter 
exit/signal port with the emerging sideband beams, and 
demodulating the resultant photodiode output signal at 
I'mod- Efficient coupling of the RF-sidebands to the 
signal port is achieved by incorporating a macroscopic 
length asymmetry between the two arms of the PRC, 
Lasymm = (^2 — L^)/2, thus allowing an optimal fraction 
of sideband power to be extracted at the exit port during 
each round-trip through the PRC, even while the carrier 
is held to a dark-fringe there. This signal generation 
method (see Sec. IIII C 4|l is referred to as the "Schnupp 
asymmetry scheme" j29| . 

To simplify our simulation task, we model only those 
aspects of the LIGO system which have a direct role in 
generating the GW-signals: the aforementioned carrier 
beam, and its RF-sidebands, resonating in the core opti- 
cal system of Fig. ^ If we make the (good) approxima- 
tion that most of the sensitivity-limiting shot noise at the 
beamsplitter exit port is due to these fields — some car- 
rier field power emerging (primarily) because of imperfect 
dark- fringe contrast, and RF-sideband field power being 
maximally channeled to the exit port — and neglect 
contributions from various other fields used for interfer- 
ometer control systems, or from other (controlled) noise 
couplings [231, then we can do a full calculation of the 
shot-noise-limited GW-sensitivity of the simulated inter- 
ferometer. Beyond this, rather than explicitly emulating 
all of the control systems of the real LIGO, the simula- 
tion program uses alternative methods (see Sec. IIII C\ for 
the numerous parameter adjustments (resonance finding. 
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etc.) that are needed to optimize interferometer perfor- 
mance; and it does so in a manner that reflects the be- 
havior of the real system as accurately as possible. 

In summary, the primary goal of our modeling efforts is 
to completely solve for the steady-state carrier and side- 
band fields that resonate in a LIGO interferometer with 
realistic optical imperfections, which is held to its proper 
operating point, and is optimally configured for maximal 
GW-sensitivity. The full effects of optical imperfections 
upon the sensitivity of interferometric detectors are then 
determined from these simulated cavity fields. 



FIG. 2: Transverse slice of a Hermite-Gaussian TEMio mode 
III. THE FUNDAMENTALS OF THE (only the real part is shown), taken at the waist plane of the 

SIMULATION PROCESS beam, and recorded on a 64 x 64 pixelized grid. 



Here we present an overview of the program and its op- 
erations during execution, and the additional tasks which 
must be done during pre- and post-processing stages in 
the course of performing simulation experiments. The 
full technical details of our work can be found in p^. 

In regards to computing platforms, the simulation pro- 
gram was originally written in the SPARCcompiler 3.0 
version of Fortran 77. A complete run (including carrier 
and sideband frequencies, with all interferometer fields 
computed and all optimizations done) with a set of non- 
ideal mirrors, and with (for example) 128 x 128 pixelized 
grid maps used for the optical computations, takes less 
than a day on a 2-processor SPARCstation 20. This com- 
putation time is reduced to a couple of minutes when the 
program is executed on massively-parallel supercomput- 
ing platforms (e.g., H^^, [13^). The runs to be pre- 
sented in this paper were performed using 32 nodes of 
the Paragon machine Trex a 512 (compute) node 
machine utilizing Intel i860 processors. 



A. Basic Optical Operations 

We use the customary approach [II [H El El IM El 

l2^ for the grid-based modeling of the laser field wave- 
fronts: a primary propagation direction is assumed for 
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UNIX cluster (skinner.cacr.caltech.edu) at the California Insti- 
tute of Technology, and for using it to conduct extensive compu- 
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form. 



the collimated beam(s) in each part of the interferometer, 
and a perpendicular slice can be taken anywhere along 
the beam propagation axis, at locations of interest. Each 
of these beam slices is recorded on a two-dimensional 
(2D) grid, with a pixel entry representing the complex 
electric field ( "e- field" ) amplitude at that transverse spa- 
tial position in the slice. Figure[5]is an example of such a 
grid-based electric field, in particular that of a Hermite- 
Gaussian TEMio mode 114]. No polarization vector is 
currently recorded in our grids (we cannot, for example, 
model birefringence effects). The precision that can be 
achieved by the program depends upon how accurately it 
simulates the two basic physical processes which must be 
performed on the interferometer e- fields: propagations^ 
and interactions with mirrors. We discuss propagations 
first. 

The program utilizes Siegman's method E3 foi' the 
plane-to-plane propagation of light in the paraxial ap- 
proximation, performed via a three step process: a 
Fourier transform and an inverse transform, sandwiched 
around a pixel-by-pixel multiplication step in spatial- 
frequency-space ("fc-space"), in which the e- field slice 
is multiplied with a distance-dependent 2D propagator 
matrix of the same pixelization. For each curved, poten- 
tially jagged mirror profile to be simulated, a flat plane 
is defined near that side of the optic to serve as a refer- 
ence plane. Propagations thus translate an e-field slice 
through the large distance from an initial reference plane 
to a destination reference plane near a mirror (or at any 
chosen position along the beam axis). The computation- 
ally intensive parts of this process (i.e., the transforms) 
can be performed rapidly by a Fast Fourier Transform 
(FFT) routine (e.g., 32]). All such macroscopic-distance 
propagations (i.e., AL ^ Acan) of c-ficlds are done in 
this manner. 

A very important issue to deal with for propagations is 
the problem of aliasing )32j |. a common complication for 
calculational procedures that use discrete Fourier trans- 
form methods. Aliasing can create artifacts in the sim- 
ulation results if (relatively) large-angle scattering sends 
power beyond the edge of the grid calculational window 
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during the propagations. Such power automatically re- 
enters the calculation from the other side of the grid, and 
may (fraudulently) be incorporated back in the physical 
simulation, instead of being filtered out, as it would be 
by absorbing baffles in a real interferometer. This prob- 
lem is most significant for the long propagations through 
the 4 km FP arm cavities. 

We can reduce (or eliminate) such aliasing by filter- 
ing the relevant fc-space matrix that will multiply an e- 
field slice, in between the Fourier transform and inverse- 
transform steps of a propagation. The goal is to preserve 
as much "real power" , while eliminating as much "aliased 
power" , as possible; the distinction between them is that 
the former always stays within the apertures of the finite- 
sized mirrors, while the latter leaves the grid completely, 
but re-enters from the other side and comes far enough 
into the middle of the grid to fall once again within the 
mirror apertures. 

Let W be the physical side-length of the square cal- 
culational window, A be the mirror aperture diameter, 
and L be the propagation distance length. The cut- 
off between physically real power and "aliasing" power 
is a matter of propagation angle: power traveling at 
6 < 0^ = A/L may be able to stay within the mir- 
rors, while power with 9 > 9a. = (W — A)/L usually 
cannot, but will sometimes be able to leave the grid 
and return to erroneously re-enter the mirror apertures. 
For a discrete grid with A'pix pixels on a side (numbered 
n — — (iVpix/2 — 1) . . . (A'^pix/2)), these angles correspond 
to (respectively) the pixel numbers TVr = Int[A-VF/(L-A)] 
and A^a = lni[{W - A) ■ W/{L ■ A)], where A is the hght 
wavelength. All of the aliasing power can be safely elim- 
inated by nulling pixels with \n\ > A^a in the fc-space 
propagator matrix, without removing any real power, as 
long as iVa > iVr. But if A^a < A'r, then one must choose 
some compromise between keeping real power and cut- 
ting out aliasing power for pixels < \n\ < N,-. One 
can, of course, force A'a > N,- to be true by increasing 
the calculational window size to be much larger than the 
mirrors (i.e., W > 2 A), but this requires more pixels to 
be used in the grid (e.g., 256 x 256), which increases the 
computational load. Making W large without using an 
adequate number of pixels would lead to poor sampling 
of the laser beam itself ,24j , thus creating a new aliasing 
problem, due to the inadequate resolution of the discrete 
grids. In our runs, we use large calculational windows 
and many pixels whenever possible (such as for the runs 
presented below); but when necessary, the program has 
an option which applies an apodization scheme to "trim" 
the propagation operators in a way that enforces a grad- 
uated compromise between keeping real power and elim- 
inating aliasing power. 

Next, we consider mirror interactions. The program 
must carry out two basic mirror interaction operations: 
reflections (r) and transmissions (t). Mirrors will not be 
perfectly uniform or flat in the direction transverse to 
the beam; they will have inhomogeneities in the refrac- 
tion index and/or thickness of their substrates, spatially- 



varying surface height proflles, variations in the quality 
of the reflective coatings, macroscopic curvatures, tilts, 
etc. This leads to optical path length variations (i.e., 
spatially- varying phase delays) across their proflles, as 
well as variations in the amplitudes of r and t, the lat- 
ter mostly due to variations in the reflective side or anti- 
reflective (AR) side coatings. These phase and amplitude 
effects are simulated by creating complex mirror maps 
for all r and t operators, which are also recorded on 2D 
pixelized grids. As demonstrated by Vinet et al. [l5l| . 
it is a good ("short distance") approximation to treat 
each pixel of the beam as an independent little plane 
wave, and reflect (or transmit) that piece of the e-fleld 
by multiplying that pixel in the e-fleld map by the corre- 
sponding pixel in the relevant mirror map, so that each 
e-field pixel interacts only with the mirror pixel located 
immediately in front of it. Thus each mirror reflection 
or transmission operation is reduced to a pixel-by-pixel 
multiplication step between a mirror operator map and 
the e-fleld slice on the reference plane near to it. The sim- 
ulation program uses these mirror interaction operations, 
in combination with the propagation algorithm described 
above, to model the entire behavior of the laser flelds in 
the interferometer. 

There are some limitations in using this methodology 
for mirror interactions, besides the obvious one of re- 
stricting it to operations over short distances (i.e., in 
the near-fleld). Vinet et al. derive a requirement 

which specifles how small the deviations of a mirror pro- 
flle can be from its idealized shape, before this pixel-by- 
pixel multiplication method loses a large amount of accu- 
racy compared to the exact calculation via the Huygens- 
Fresnel integral formulation in scalar wave theory |14| . 
Realistically imperfect (LIGO-quality) mirrors easily sat- 
isfy this requirement. Furthermore, as noted by Tridgell 
et al. |l7j |. the difference in phase between one pixel in a 
mirror map and its neighbor must be smaller than 27r if a 
continuously- varying surface displacement on the mirror 
is to be adequately sampled by the grid. For our sim- 
ulation parameters (cf. Table P) , this limits mirror tilts 
io 9 < 10"'' radians (not including the base 45° tilt of 
the beamsplitter, which is handled separately), and mir- 
ror curvature radii to -Rcurv > 1-3 km; both of these 
limitations are easily satisfied in our runs. Finally, we 
supplement these requirements with a general rule-of- 
thumb: a tiny Gaussian beam, with a waist size equal 
to the width of a pixel [and an initial propagation di- 
rection with respect to the beam axis that is defined by 
the e-field's wavefront curvature at the location of that 
pixel), must not expand or shift over so much so that it 
would relocate a significant fraction of its power onto any 
neighboring pixels, during the entire course of the mir- 
ror interaction. If that rule is broken, then this method 
for mirror interactions will be inaccurate, and the pixcl- 
by-pixel multiplication of maps will not be sufficient for 
modeling reflection or transmission operations. 

To create a complete (but not over-determined) de- 
scription of a mirror, we consider the full information 



6 



expressed by mirror interactions, as well as physical sym- 
metries and the conservation of energy. A 2-port mirror 
(reflective- and AR-sides) must relate 2 complex input 
fields to 2 complex output fields. Thus 4 complex (or 
8 real) elements are needed, at each pixel location, to 
specify that pixel of mirror completely. These 4 com- 
plex elements correspond to the complex reflection and 
transmission operations performed from the two differ- 
ent sides of the mirror. The two transmission operations 
from either side — barring excessive beam expansion or 
focusing during the transmission — must in fact be the 
same |33|: tright = iieft = t (except for the beamsplitter, 
which requires distinct transmission maps for the inline 
and offline paths). The mirror description is thus reduced 
to 3 complex elements per pixel, i.e., 3 complex mirror 



P before nafter a | jribefore \2 , a | jribefore |2 

- = Aj-cR. ■ lii'rofi.-sidcl + ^AR ■ I^AR-sidcl 



For the simplest case of a loss-free mirror (^rofl. = ^ar 
= 0, and pbeforc _ paftcr ^ ^^lis just rcduccs to the 
complex generalization of one of Stokes relations [33l |: 
{rrcfi. / r'^B.) = ~(*/^*)- Assuming the transmission coef- 
ficient t to be real, this generates the familiar possibil- 
ities for the phase relationships between the reflectivi- 
ties on either side of a mirror: e.g., (rrofl.,7'AR) oc 
(?'rcfi., ^ar) oc (±1, =f1), etc. For the more general case of 
a lossy mirror, we can convert Eq. |^ into a simpler pre- 
scription by considering variations to the relative phases 
and amplitudes of E^^l°J^^^^ and -Ear^mc ^^^^ generating 
the strict energy conservation condition: 

, ^-^'-"ft-'-^AR ^2) 

\t* ■r,,fi,+t-rl^\ - ^ ' 

This inequality (generalized further for the beamsplitter) 
must be satisfied (at every pixel) to guarantee conserva- 
tion of energy in a mirror, given any e-fields which could 
be incident upon it. Our simulation program enforces 
this condition by testing the input data maps for each 
mirror, and rejecting runs with unphysical specifications. 

Given all of these requirements, the pixelized mirror 
maps described above are versatile tools for modeling a 
wide variety of optical features and imperfections, as dis- 
cussed in Sec. 

B. Relaxation of Steady-State Electric Fields 

Our program models a static interferometer, assumed 
to be held to the correct operating point for an indefinite 
period of time. This assumption enormously simplifies 
the program, and our overall simulation task, while still 
enabling us (see Sec. IIII D|) to compute the frequency- 
dependent GW-response of a LIGO interferometer. (The 



maps: the transmission map and reflection maps from 
either side, rj-cfl. and tar. 

Next, by energy conservation we have (for each mirror 
pixel): l-|i|2-|n.efl.|2 = > 0, and 1- |rARp = 

^AR > 0, where ^refi., ^ar are the losses experienced in 
reflection from either side of the mirror. But these con- 
ditions are not sufflcient. By considering the complete, 
superposed e-fields which exist on either side of the mir- 
ror, both before and after the mirror interaction (i.e., the 
incoming e-fields vs. the outgoing ones), and by requir- 
ing that the total power in all e-fields must not increase 
due to the interaction, one may obtain the following in- 
equality (expressed in terms of the incoming, "before" 
fields): 



2g?{(r . r,efl. + t ■ rlR,) . E^±!l^, ■ (i^^f ^f^c)*} > • (1) 
I 

dynamics of LIGO control systems do not necessarily 
come into this calculation, since the mirrors act like 
"free-masses" at the GW-frequencies that LIGO is most 
sensitive to, and because the GW's are presumably not 
strong enough to interfere with the interferometer's reso- 
nant lock.) The primary task of the simulation program, 
therefore, is to compute the relaxed, steady-state, reso- 
nant electric fields that build up in each part of the in- 
terferometer, when it is excited by a laser beam of fixed 
amplitude, orientation, and frequency (at the nominal 
carrier and sideband frequencies), entering through the 
power recycling mirror. 

The multiply-coupled-cavity nature of the interferome- 
ter, and the physical complexity of the simulation, means 
that the program must iteratively solve for an irreducible 
number of resonant e-fields (one per cavity) in the inter- 
ferometer. For the first-generation LIGO configuration, 
three e-fields must be computed via relaxation: one in the 
PRC, and one in each of the FP arm cavities (though the 
precise location of each relaxed e-field in its cavity can 
be arbitrarily chosen) . When the relaxation algorithm is 
finished, these three principal e-fields can be propagated, 
refiected, transmitted, and/or superposed together in or- 
der to generate the complete steady-state e-fields every- 
where in the system. 

For each e-field that will be relaxed, there is a steady- 
state equation of the form: 

-^steady-state — A^E/gtcady-state } -^exc ■ (3) 

The left hand side of Eq. Q is the e-field to be solved for; 
we display it here as a vector because it has a propaga- 
tion direction: either "forward" or "backward" along the 
beam axis. The expression, A{i?steady-state}, represents 
this steady-state e-field after it has gone on a round-trip 
through the interferometer, and has returned to its start- 
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ing point. (Note that this "round-trip" includes all possi- 
ble closed-loop paths through the interferometer that do 
not ever pass through the same location — with the same 
propagation direction — of any of the other principal e- 
fields being relaxed.) Lastly, iSexc represents a composite 
excitation e-field, which may consist of both a primary 
excitation e-field (e.g., the input laser beam, for the e- 
field being relaxed in the PRC), plus any "leak" e-fields 
that arrive from the other principal fields being relaxed 
(e.g., fields leaking from the FP-arms into the PRC). It 
is this leak-field part of the excitation term which is the 
source of many coupled-cavity effects in the interferom- 
eter; and it is the presence of the A{i<^stcady-statc} term 
in Eq. (0) which makes i?stcady-statc "self-coupled" , thus 
requiring us |3| to use an iterative process to solve for 
it. 

Specifying a relaxation convergence scheme is done by 
prescribing how the {N + 1)**^ iterative guess for a given 
e-field is obtained from its iV**^ iteration (and from the 
N^^ iteration of the other e-fields being relaxed simulta- 
neously). Considering Eq. (|3J|, the simplest possibility is 
to make the choice: 

B^+i) = A{SW}+iW . (4) 

This process can then be repeated for as many iterations 
as necessary, until the steady-state equation, Eq. 

is satisfied by i?'^) to within a pre-specified threshold 
of accuracy (typically 1 part in 10^ in power, for our 
runs). This relaxation formula is guaranteed to succeed 
at smoothly converging an e-field to its correct steady- 
state form — barring complications which may arise 
from changes to the interferometer caused by the param- 
eter optimization routines (see Sec. IIII C|l — because 
it imitates how power actually does build up, through 
a sum of many bounces, in the cavity system of a real 
interferometer. This iteration process has been used suc- 
cessfully by previous researchers (e.g., fisfl. 

Since this method does model the true physical buildup 
of power, however, it requires a great many iterations 
to converge, especially for coupled-cavity systems with 
large Q-factors (i.e., long transient decay times). We have 
found this relaxation scheme to be forbiddingly slow for 
a fuU-LIGO simulation program. 

Therefore, as first suggested (and implemented, for a 
simpler cavity arrangement) by a LIGO colleague [33| . 
our simulation program uses a different approach. In- 
stead of Eq. ^ for choosing the {N + 1)*"^ iteration, we 
use the following expression: 

^a.iW+6.(A{iW}+iW) + ,.iW , (5) 

where (a, 6, c) are unknown, complex coefficients that are 
solved for by minimizing the error in what the steady- 
state equation will become in the next round, with 
£■(^+1) (as a function of these unknown coefficients) tak- 
ing the place of i^stoady-state in Eq- ®, and E^^c taking 
the place of Ecxc (noting that Eixc does not yet take into 



account the new (a, b, c) coefficients that will be chosen 
for the other e-fields being iterated). The resulting ex- 
pression for the steady-state error can be differentiated 
with respect to the six available degrees of freedom (the 
real and imaginary parts of a, 6, and c), resulting in six 
simultaneous equations that are solved via matrix inver- 
sion. This relaxation method (which essentially reverts 
to the simpler method if we set (a, b, c) = (0, 1, 0)) gains 
a huge advantage compared to the method of Eq. Q , by 
having these additional, useful degrees of freedom avail- 
able for each iteration [SJI ■ The number of iterations nec- 
essary to achieve convergence are greatly reduced (often 
by ^^1-2 orders of magnitude), resulting in much faster 
e-field relaxation. 

This "abc" iteration algorithm does have certain draw- 
backs. It is more difficult to implement in code, and it 
requires more information to be stored (thus using more 
memory), and more propagations to be performed (de- 
spite what is stated in 34], because we use three un- 
known coefficients instead of two: c 7^ 0), during each 
iteration. It is also somewhat less stable, because of the 
wide range of abc-coefficients that may be chosen by the 
error minimization algorithm; in fact, large excursions 
in the iterated e-field structure and power level appear 
to be required during (typically) the first ~50-100 itera- 
tions, in order for the accelerated convergence scheme to 
function properly. We have found that the best (straight- 
forward) way to enhance the convergence stability of the 
"abc" relaxation algorithm, without slowing it down to 
the pace of the method of Eq. Q), is to place hard limits 
upon the allowed choices of the abc-coefficients during 
the early stages (the first ^100-200 iterations) of runs. 
This allows the relaxation process to avoid failure during 
the initial, large e-field adjustments, after which it settles 
down to converge efficiently to the steady-state solution. 

Lastly, we note that the choice of (a,6, c) for any one 
field will affect future iterations of the other fields, be- 
cause they are all coupled via the Ecxc leak-fields. Iter- 
ating them independently from one another ignores this 
coupling, and causes some slowing (and oscillation) of 
the relaxation process (we see a temporary "sloshing" of 
power back and forth between the PRC field and the FP 
arm cavity fields). This problem can be eliminated by 
choosing the (a, b, c) coefficients for all three relaxed e- 
fields simultaneously, by calculating the 18 unknown pa- 
rameters (i.e., 9 complex "abc" coefficients) in order to 
minimize a specially- weighted sum of the iteration errors 
for all three relaxed fields. Such a "global abc" relax- 
ation scheme is significantly more demanding in terms 
of coding difficulty and computer memory usage, but it 
is capable of making the relaxation process more sta- 
ble, and further reduces the number of iterations needed 
to reach convergence. Though we have not yet imple- 
mented this global relaxation scheme into simulations of 
the first-generation LIGO interferometer, we have suc- 
cessfully incorporated global abc relaxation (with four 
relaxed e-fields needed, instead of three) into the version 
of our code used for the advanced-LIGO Dual Recycling 
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configuration, with good results |24||. 

In summary, the "abe" iteration method presented 
here is a sufficiently reliable and extremely fast relax- 
ation scheme which greatly reduces program execution 
time, thus enabling us to perform fuU-LIGO simulation 
runs (with significant optical deformations) in a reason- 
able amount of time. 



C. Parameter Optimizations 

At the core of our efforts to make a realistic simu- 
lation of a LIGO interferometer are several procedures 
which bring the system into an "optimal" configuration 
for signal detection. The problem of optimization is a 
highly nontrivial matter: not only does the interferom- 
eter possess numerous degrees of freedom that must be 
optimized (e.g., all resonant cavity lengths, the Schnupp 
length asymmetry, etc.), but each evaluation of perfor- 
mance (i.e., GW-sensitivity) as a function of these opti- 
mizable parameters is extremely time consuming, since it 
requires the detailed computation of the carrier and side- 
band e- fields everywhere in the interferometer. It there- 
fore appears infeasible to use the brute-force method of 
optimizing the interferometer's GW-sensitivity function 
by evaluating it for a thorough sampling of points over 
the entire, multidimensional parameter space. 

As an alternative, each key parameter is optimized 
separately in our program, using some error signal (to 
be brought to zero), or some function of merit (to op- 
timize), which is strongly (and solely) dependent upon 
that individual parameter, and which is a true measure 
of when that parameter is well chosen for the maximiza- 
tion of GW-sensitivity. This strategy is aided by the fact 
that the interferometer's GW-sensitivity is quite insen- 
sitive to particular combinations of parameter changes, 
such that if certain optimizable parameters are displaced 
from one apparently optimal point in the multidimen- 
sional parameter space, then the other parameters will 



adjust themselves (via our optimization procedures) to 
compensate with virtually no reduction in overall inter- 
ferometer sensitivity. Wc note that these optimization 
routines in our program are performed concurrently with 
the e-ficld relaxations, so that the final results emerging 
from the iterative relaxation scheme are the steady-state 
e-fields of a fully- optimized interferometer. 

The following subsections give an overview of the pa- 
rameters that the program optimizes, the criteria for op- 
timizing them, and the physical considerations which un- 
derlie their significance. Note that we have performed 
extensive empirical tests to verify that each of the opti- 
mization procedures discussed below works as desired. 



1. Length Adjustments for Carrier Resonance 

As discussed in Sec. ^ the carrier frequency beam 
must have a double resonance in the system consisting of 
the power recycling cavity and Fabry-Perot arm cavities, 
while being held to a dark-fringe at the exit port of the 
beamsplitter. Our method for achieving these resonance 
(and dark-fringe) conditions is similar to that of Vinet et 
al. 15], in which we null the computed "phases" between 
certain specified cavity e-fields by adjusting the various 
cavity lengths. (A length change will alter the phase re- 
lationship between an e-field that has taken a round-trip 
through an adjusted path length, and one that has not 
— potentially bringing a cavity to resonance.) Alterna- 
tively, we have chosen not to use the method of McClel- 
land et al. in which the e-fields are re-relaxed for 
each trial set of lengths until the configuration for max- 
imum power buildup is found, since that would involve 
a time-consuming search over a multi-dimensional phase 
space of independent, controllable cavity lengths. 

The phase $ between two e-fields is defined via an 
overlap integral (or more precisely, by a discrete sum over 
the pixelized N x N grids), as follows: 



^[g„g.].tan-M ^<g^'g-^ ], 
< Ei\E2 > 

(Gale. Window size) ^ t^^, s , s 

with: < E,\E2 >= ^ ^ • E ^^(*'^) • ^2(^, j) ■ (6) 

i=l 3 = 1 



The phase between an e-field and its round-tripped ana- 
logue can be driven to zero by a (sub-wavelength) path 
length change, computed via the formula: AL — — $ • 
(— Aiascr/47r) = $/2fc (an extra factor of 1/2 is included 
here because the round-trip doubles the phase change 
caused by a cavity length adjustment). Note that this 
procedure cannot simply be performed once — these 



phases depend in detail upon the structure of the iter- 
ated e-fields, and they must be repeatedly measured and 
adjusted throughout the e-field relaxation process. 

To achieve resonance in a particular cavity, the 
proper solution is to ensure that the fields £'stcady-statc , 
A{£'stcady-statc}, and -Eoxc (as defined in Sec. IIIIB|) all 
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have zero phase between them. (This condition of mu- 
tually zero phase is g iven some theoretical justification 
in the literature |34j . as well as being apparent from 
Eq. ©)• The phase- nulling procedure is performed in 
the FP arm cavities and in the PRC by performing mi- 
croscopic adjustments to the three relevant cavity lengths 
here: L4, L5, and the "common-mode" recycling cav- 
ity length, LpRc- Similarly, the dark- fringe condition 
at the beamsplitter exit port is achieved by setting the 
phase between the carrier e-fields coming from the in- 
line/offline recycling cavity arms to an odd multiple of 
TT, by microscopically adjusting Lasymm = {L2 ~ Lj,)l2 
via "differential-mode" corrections to L2 and L^. 

Unlike the procedure discussed in we do not 

choose any particular spatial mode (such as the lowest 
order, TEMqo mode) of the interferometer e-fields for cal- 
culating these phases. Rather, we use the entire e-fields, 
for two reasons: first, for the dark-fringe condition, one 
wishes to minimize the total dark-port power emerging 
from the beamsplitter exit port, all of which contributes 
to the shot noise; second, for the carrier resonance con- 
ditions, coupling between modes brings power back into 
the TEMqo mode from higher modes, so that bringing 
the total field (i.e., the "perturbed interferometer mode") 
to resonance also maximizes the TEMqo power used in 
calculations of the GW-signal. In any case, we have ob- 
served that these two methods (i.e., with or without spa- 
tial mode selection for resonance-finding) typically give 
very similar results. 



2. Sideband Frequency Fme-Tuning 

The RF-sidebands must also satisfy important condi- 
tions, specifically resonance in the PRC, and near-anti- 
resonance in the FP arm cavities. These conditions are 
affected by the optics, and the sidebands must also there- 
fore be tuned for optimum performance. 

As a first approximation, the cavity lengths and RF- 
modulation frequency are initialized through an analytic 
evaluation designed to simultaneously optimize carrier 
and RF-sideband performance. But finer adjustments 
are required for the PRC resonance condition, and since 
the cavity lengths are already fixed by the resonance re- 
quirements of the carrier beam (the carrier and sideband 
e-fields are relaxed in separate code executions, with the 
carrier first), the free parameter which remains to be ad- 
justed is the sideband frequency. In a manner analogous 
to that specified for cavity length changes, sideband fre- 
quency adjustments for resonance are periodically com- 
puted as: Aj/sb = ^ ■ (c/47ripR,c); though somewhat 
smaller frequency changes are actually made in practice, 
during each adjustment, to lessen the disturbances to 
the e-field relaxation process. The cumulative frequency 
change for typical runs is small (usually a few hundred 
Hz or less), but necessary to achieve sideband PRC res- 
onance. Lastly, we note that the overall results of runs 
for the two different RF-sidebands (i.e., the "upper" and 



"lower" sidebands, ±vsb = ^'carri^'mod) are usually fairly 
similar; we generally perform computations for only one 
sideband, and assume mirror-image results for the other. 



3. Recycling Mirror Reflectivity Optimization 

The level of power buildup in the interferometer, and 
hence its GW-response, depends upon the reflectivities 
of the mirrors, and is constrained by the mirror losses. 
In LIGO, the reflectivities of most of the mirrors have 
been predetermined according a variety of auxiliary phys- 
ical requirements (e.g., p^). But the principal criterion 
for specifying the reflectivity of the power recycling mir- 
ror (i?i) is the maximization of the gain that is achiev- 
able with power recycling. The power recycling gain, in 
turn, depends critically upon the losses experienced in 
the imperfect interferometer, which cannot be precisely 
determined via analytical estimates. Ri is therefore a 
parameter that can be optimized by the program dur- 
ing each run'^, in order to determine the best achievable 
recycling gain, given the particular optical imperfections 
being studied. The reflectivity of a real mirror, of course, 
is not something that can ordinarily be adjusted on the 
fly; rather, the results of this optimization routine in our 
simulation program have been used to help determine 
appropriate "design values" of for power recycling 
mirrors that are procured by LIGO. 

For interferometer losses that are small, it can be 
shown (e.g., [s^) that the optimal choice for the recy- 
cling mirror transmission is Ti^optim ~ ^ifOj where ^ifo 
is the effective total loss in the full interferometer (includ- 
ing the loss in the power recycling mirror). The optimal 
reflectivity is therefore given by: i?i,optim = l-Ti.optim ~ 
1 — (^1 — AiFo)- The value of ^ifo, and thus of i?i,optim, 
will be strongly affected by optical deformations. 

The quantity used for this optimization procedure is 
the real part of the integrated overlap between the imme- 
diate ( "prompt" ) reflection of the laser excitation beam 
from the AR-side of the power recycling mirror, with the 
complete, composite e-field that is ultimately reflected 
back from the interferometer. (The imaginary part of 
this overlap integral is related to interferometer reso- 
nance, and will essentially be zeroed if the cavity lengths 
are properly adjusted for carrier resonance in the interfer- 
ometer, as per Sec. IIII (TT\ ) The real part of the overlap 
constitutes an error signal^ that can be driven to zero 



Rl optimization is performed specifically during carrier runs, 
since it is the carrier which requires a high PRC gain; the side- 
bands do not benefit from long PRC storage times, and would 
in fact suffer less degradation with a PRC storage time of zero, 
and immediate ejection at the beamsplitter exit port. 
* This error signal is only precisely valid when the power re- 
cycling mirror obeys the "lossless-mirror" Stokes condition, 
(''refl. Aar) ~ i^^- Sec. IIII A1 . Our program automat- 

ically obeys the amplitude part of the condition (i.e., |rj(,fl.| = 
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by changes to Ri. The magnitude of each change in the 
optimization process will be proportional to this error 
signal, though the proportionality constant (which we 
have selected through empirical tests) is not crucial, as 
long as it is large enough to achieve rapid optimization 
(and close convergence to the optimal value) , while being 
small enough to ensure a stable optimization process. 

This nulling of < £^prompt|^^totai reflection > IS equiv- 
alent to minimizing the total ("interferometer mode") 
power that is reflected from the system, thus maximizing 
the power that is dissipated (and hence, that is circulat- 
ing) inside of it. In practice, LIGO has chosen an initial 
value for i?i that is slightly below such an "optimized" 
value, both to hedge (on the safer side) against uncer- 
tainties (and gradual increases with time) of the effective 
interferometer losses, and to provide some reflected light 
for length control signals. For the specific runs to be pre- 
sented in this paper, the recycling mirror reflectivity has 
been driven all the way to i?i,optim; though runs can also 
be performed with our program in which i?i is held to 
any particular fixed value that may be preferred due to 
practical considerations. 

4- Schnupp Length Asymmetry Optimization 

Also noted in Sec.HTIwas the incorporation of a macro- 
scopic length asymmetry (~few tens of cm) between the 
inline and offline paths of the PRC, in order to max- 
imally channel sideband power out through the beam- 
splitter exit port, for use as a local oscillator in the het- 
erodyne GW-signal detection scheme. The maximization 
of this local oscillator light requires a careful balance be- 
tween extracting sideband light from the interferometer 
promptly, before signiflcant power is wasted due to mir- 
ror losses; yet leaving the sidebands in the interferometer 
long enough (i.e., for enough round-trip bounces) so that 
they can take full advantage of the broadband amplifica- 
tion provided by power recycling. 

The total phase between the inline and offline RF- 
sideband fields, when they meet at the beamsplitter, is 
analytically given as: $asymm = 2 x [-2 k^od iasymm], 

I^'arI) recycling mirror; and though the phase part of 

it is not mandated, we impose it for virtually all of our runs, 
including those described below in Sec. IIVI 
® Eq. J7J is actually a simplification of the proper formula; the 
simulation program uses a more complete expression |24| , which 
does not assume balance between the two paths or a 50%-50% 
beamsplitter, and which the program must solve numerically. 



where fcrnnH and Xasymm are defined as in Sec.lHJ It can 
be shown j2J| that the maximum local oscillator light is 
generated by the choice: 

COs(<I>asyinm) = 2^i?i i?arm avg. ^BS , (7) 

where i?arm avg. is the average refiectivity experienced by 
the beam along its inline and offline paths (including all 
Fabry-Perot arm cavity effects), and where Rbs ~ ^bs 
for the beamsplitter has been assumed here^. 

The reflectivity values to be used in Eq. ||7|) (or in 
its generalized version for an unbalanced interferometer) 
cannot be analytically determined before the program 
is run, since i?arm avg. and Rbs depend upon the losses 
due to optical imperfections, and i?i is a variable that 
is optimized during the carrier run. The simulation pro- 
gram therefore uses the resonating powers at various in- 
terferometer locations (during the RF-sideband run) to 
determine "effective" reflectivities for use in Eq. JTI), to 
estimate the optimal choice of ^' asymm- Periodic changes 
to iasymm (totaling a few cm, cumulatively) are per- 
formed during the sideband field relaxation iterations, 
until $asymm IS brought to this desired value. These 
macroscopic length changes are applied antisymmetri- 
cally to the inline and offline paths, and in integral mul- 
tiples of Acarr/4, in Order to avoid (as much as possi- 
ble) any disruption to the carrier resonance/dark- fringe 
conditions®. Lastly, we note that the program considers 
only the TEMqo components of the e-fields while eval- 
uating <I>asymm and Eq. JTJ) (as opposed to considering 
the total e-fields for optimization, cf. Sec. IIII (JT]) . since 
this extracted RF-sideband light is to be used directly 
for generating the GW-signal. 



5. Sideband Modulation Depth Optimization 

The sideband fields are generated via electro-optic, ra- 
dio frequency modulation of the initial, carrier-frequency 
laser field. This process is represented mathematically as 
follows: 

The program also uses measured values of "3? asymm, not this an- 
alytical one, in its Lasymm-optimization adjustments. 
® We typically follow the initial carrier/sideband 2-run set with an 
additional carrier/sideband 2-run set, with no changes to Lasymm 
in the latter set, to verify that the proper conditions for the 
carrier have not been disturbed. 



oo 

Mji/ias e } = Ei^s e ^ = i/ias e • J„(r) e 

71— — OO 

« i;ias • {Jo(r) e*"*-f Ji(r) p5'"+^)* + j_i(r) e^f"-^)*} , (8) 

I 

where: i?ias is the (constant) laser field amplitude, lu = the modulation angular frequency, F is the modulation 
27ri'carr IS the carrier angular frequency, S = 27ri^niod is 
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depth, J„(r) = (— l)"J_„(r) is the Bessel function (with 
integer order n) as a function of F, and where we have 
dropped the higher-order terms in the series because T 
wiU be kept smah enough for them to be unimportant. 

Specifying the modulation depth is a matter of dividing 
up the total available laser power between the carrier and 
its RF-sidebands. The carrier and sidebands both con- 
tribute to the GW-signal, and due to their (unavoidable 
for the carrier, and deliberate for the sideband) ejection 
at the beamsplitter exit port, they both contribute to 
the shot noise that competes with this signal. The opti- 
mal modulation depth is determined via maximization of 
the overall ratio of GW-signal to shot noise. This ratio 
(given explicitly in Sec. IIII Dl below) is dependent upon 
the buildup of the steady-state e-fields, and hence upon 
the specific interferometer deformations which exist; for 
a real LIGO detector, in particular, the best modulation 
depth (or the "best achievable" F, given technical limi- 
tations) to impose on the carrier laser beam can only be 
precisely determined once the differing performances of 
the carrier and sideband fields in the interferometer are 
known. 

In our modeling work, the modulation depth opti- 
mization is performed as a post-processing step, after 
the simulation runs for the carrier and sideband e-fields 
have been completed. Besides the straightforward op- 
timization of the shot-noise-limited GW-sensitivity (see 
Eq. Hll|) below) with respect to F, we are also careful to 
make sure of two things: that F is small enough so that 
the higher-order modulation terms are indeed unimpor- 
tant as far as their potential buildup in the LIGO interfer- 
ometer is concerned; and that the total amount of power 
(carrier plus sidebands) falling upon the GW-signal de- 
tection photodiode at the beamsplitter exit port is not 
excessively large. 



6. Mirror Tilt Removal 

The final optimization procedure discussed here is a 
pre-processing step that is performed upon the mirror 
profile maps before they are read in during execution of 
the simulation program. Specifically, it is the removal 
of any overall tilts due to mirror surface variations or 
substrate inhomogeneities. (A "substrate tilt" is merely 
a mirror thickness wedge, which can be countered in a 
real interferometer by a small change to the reference 
axis of a cavity. ) 

For irregular mirror profiles, which do not possess a 
unique definition of tilt, we choose the most useful def- 
inition for a LIGO interferometer: the tilt that is expe- 
rienced (or "weighted") by a Gaussian-profile incident 
beam. To first order, the effect of mirror tilts about 
the two axes perpendicular to a (predominantly Hermite- 
Gaussian TEMqo) incident beam is to generate power 
in the TEMin and TEMqi modes, in a reflection from 
that mirror j^]- These modes will have imaginary am- 
plitudes (assuming a real incident beam amplitude) if 
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FIG. 3; Sample map of a mirror surface with realistic defor- 
mations, after being processed by the tilt-removal algorithm. 
The mirror map, with most of the border region (lying be- 
yond the finite-sized mirror apertures) clipped here for visual 
clarity, is oriented such that the incident beam propagation 
axis is along the vertical direction. The width of the grid 
shown here is ~ 25 cm (deformation heights not shown to 
scale vs. width). 



the mirror has no overall "piston" displacement. Given 
a mirror surface/substrate "height" function Z{x,y), 
and defining M{x,y) = exp[—2ikZ{x,y)] (exponentiat- 
ing each pixel individually, not the whole matrix), we 
can remove the beam-weighted tilts from a given mir- 
ror deformation map by performing small angular cor- 
rections that set 3 < TEMio|M(a;, ?/)|TEMoo >= 3 < 
TEMoi|M(a;,?/)|TEMoo >= 0. This is done after first 
removing any piston offset from the mirror, via uniform 
displacements that set 3 < TEMoo|M(x, ?/)|TEMoo > to 
zero. 

An important point about this procedure is that the 
appropriate beam spot size must be used for the modes 
in the overlap coefficients given above, in order for the 
"beam-weighting" of each mirror's tilt to be correct; but 
since the beam spot size is different at different locations 
in the interferometer, one must therefore know which spe- 
cific mirror a given deformation map will be used for, be- 
fore its tilt can be properly removed. Also, we note that 
this tilt-removal process will actually give the full mirror 
a nonzero tilt, overall; but the center of the mirror (i.e., 
the part most sampled by the beam) will be essentially 
flat (on average) with respect to the plane transverse to 
the beam propagation axis. An example of such a tilt- 
removed surface deformation map is shown in Fig. |3| 



D. Program Output Quantities and the 
GW-Sensitivity Function 

A full set of program runs results in a detailed specifi- 
cation of the final, steady state of the interferometer. Of 
the large amount of output data available (either directly 
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or after some post-processing calculations) , here are some 
of the most important quantities that we examine: 

• The relaxed powers (for the carrier and RF-sideband 
fields) at several key locations in the interferometer. 

• The complete, relaxed e-fields at all desired loca- 
tions — available for graphical display, and/or for com- 
plete modal decomposition into Hermite-Gaussian TEM 
modes, which assists in the interpretation of interfer- 
ometer conditions (e.g., TEMiq/TEMqi power indicates 
residual tilts, TEM20/TEM02 power indicates beam mis- 
match, etc.). 

• The carrier contrast defect, which quantifies how 
well a carrier dark-fringe was achieved at the beamsplit- 
ter exit port for the imperfect interferometer. Defining 
the carrier powers emerging from the relevant beamsplit- 
ter ports as Pbright and Pdark, the contrast defect is given 
as: 



1 - Contrast = 1 - C = 1 



ht 



ark 



ht 



dark 



bright 



. . (9) 

A large contrast defect implies a substantial carrier power 
loss at the beamsplitter (and thus a broadband power loss 
in the system), as well as a large carrier contribution to 
the shot noise at the signal port photodetector, and an 
excess of raw power falling on that photodetector. 

• The optimized interferometer parameters, as de- 
scribed above in Sec. IIII CI In addition to their role in 
optimizing the performance of the simulated interferom- 
eter, the computed values of these parameters often have 
intrinsic importance in terms of LIGO design considera- 
tions. 

• The GW- strain- equivalent shot noise spectral den- 
sity, /isn(/), of a single LIGO interferometer. This cru- 
cial output function allows us to directly evaluate the 
sensitivity of the LIGO detector to astrophysical sources 
of gravitational waves, given interferometers with realis- 
tic optical deformations (and a realistic heterodyne GW- 
detection scheme). 

Consider Fig. 01 which shows the three most signifi- 
cant noise sources for the first-generation LIGO interfer- 
ometers [2^: seismic, thermal, and shot noise. They 
are plotted versus GW-frequency /, as spectral densi- 
ties expressed in terms of the gravitational-wave Fourier 
amplitudes, h(f), that would induce equivalent signals 
in a LIGO interferometer. These particular curves are 
the first-generation LIGO requirements 36] for the max- 
imum contributions that would be acceptable from each 
of these three main noise categories^. The total noise 
envelope can be obtained from these individual curves 
by adding them together in quadrature (i.e., incoherent 
addition of uncorrelated noise is assumed). 



Note that the thermal noise curve in Fig. [factually represents an 
approximate conglomeration of mirror internal vibration noise, 
suspension pendulum noise, and other technical noise sources. 
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FIG. 4: Requirement curves for the primary noise sources 
expected to limit tlie GW-sensitivity of first-generation LIGO 
interferometers . 



Of these contributions, seismic and thermal noise are 
random forces which will push the LIGO mirrors around 
in imitation of GW's. Photon shot noise, on the other 
hand, is a form of sensing noise, representing the quan- 
tum mechanical limit on the accuracy to which the mir- 
ror positions can be measured, given the finite amount 
of carrier power resonating in the FP arm cavities, and 
the finite amount of local oscillator sideband light avail- 
able at the signal port. While the quality of interfer- 
ometer optics has little effect upon the level of ran- 
dom force noise contributions (other than radiation pres- 
sure noise, which should be unimportant for the first- 
generation LIGO interferometers |2a|)j the quality of the 
optics does have a direct impact upon the sensing noise 
limitation to LIGO's GW-sensitivity. The presence of im- 
perfect optics not only reduces the amount of circulating 
interferometer power available for sensing mirror posi- 
tions (and thus GW- induced mirror motions), but also in- 
creases the amount of unwanted power at the beamsplit- 
ter exit port — such as carrier contrast defect power, and 
RF-sideband power in non-TEMoo modes — which con- 
tribute to the shot noise, but not to the GW-signal. We 
therefore focus upon the shot-noise-limited region of the 
LIGO noise envelope in evaluating the effects of optical 
imperfections. For each set of output results, /isn(/) can 
be computed, and compared either to the first-generation 
LIGO requirements (see Sec. IIVI below') . or to astrophys- 
ical predictions (see Sec. 1^, in order to determine the 
effects of optical deformations upon LIGO's ability to 
detect GW's of reasonable, anticipated strengths. 

A full derivation of the formula for /isn(/) is given 
elsewhere |23|- Here we present the resulting expression, 
in terms of the relevant output data from our simulation 
program. The shot noise sensitivity limit is expressed in 
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terms of the strength of the GW needed to produce a 
signal that could match this noise. For a monochromatic 
gravitational wave of the form @ : 

h^'^if) = %/2 /icos(27r/i + 0o) , (10) 
which is incident upon a single interferometric detector 
I 



where: Pias is the total excitation laser power (before ra- 
dio frequency modulation) in Watts, hp\ is Planck's con- 
stant, z/carr is the Carrier frequency, rj is the quantum effi- 
ciency of the photodetector at the signal port, Jo(r) and 
Ji(r) represent the division of laser power between the 
carrier and either one of its RF-sidebands (cf. Eq. JHJ), 
and terms like "^'carr, exit port" ' ^^c, represent the di- 
mensionless relaxed power value (in either the TEMqo 
mode or the total in all modes, as indicated) that is re- 
ported by the numerical simulation code for a carrier or 
sideband e-field in the indicated interferometer location. 
These "dimensionless" power values for the simulated e- 
flelds are all normalized in the program to an input car- 
rier/sideband laser beam power of 1 Watt. Note that we 
have not included both the upper and lower sidebands 
separately here; in many cases it is sufficient to plug 
the simulation results for either of them into Eq. (|ll|l . 
and assume the interferometer performance to be fairly 

I 



where L is the length of the FP arm cavities, c is the 
speed of light, and a properly- weighted average has been 
performed over the results for the inline and offline FP- 
arms. We note that the dependence of /isn(/) upon 
the GW-frequency is contained within the expression 
•\/l -I- (47rTs/)^, so that it has two basic regimes sepa- 
rated by the "knee" or pole frequency, /poie = l/(47rrs), 
as follows: 

^sn(/)|/«/po,„ ~ constant , /isn(/)|/»/p„i<, oc / . (13) 



with optimal incidence angle and polarization, a GW- 
amplitude oi h = h^^^f) would produce a unity signal- 
to-shot-noise ratio when sampled over unity bandwidth 
(i.e., 1 second integration time). With these definitions, 
we have: 



(11) 



symmetrical about the carrier frequency for these RF- 
modulated fields. 

Some remaining quantities in Eq. (|11|) (e.g., 
''FP arm back, ^Bs) SiiB estimated mirror (amplitude) 
refiection and transmission coefficients for the path 
that any GW-induced signal fields would take from the 
FP-arms (where they are physically generated) through 
the interferometer, until they exit at the beamsplitter 
signal port. Finally, the quantity Tg is the effective 
storage time of the GW-induced signal fields in the 
realistically-deformed FP arm cavities. Since the explicit 
simulation of the buildup of these signal e-fields in 
the arm cavities would involve an additional set of 
e-field relaxation procedures that are not performed 
by the program, the effective storage time of these 
GW-induced e-fields in the (imperfect) arm cavities 
must be approximated, as follows: 



(12) 

I 

These two regimes are evident in the plot of /isn(/) that 
is included in Fig. ^ 



IV. A SELECTION OF RESULTS OBTAINED 
WITH THE SIMULATION PROGRAM 

In this section, we will demonstrate the code with a set 
of runs incorporating measurement maps made from very 
high quality mirrors. The only interferometer imperfec- 
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tions included are these mirror deformation maps, the fi- 
nite sizes of the mirrors, and a small amount of pure loss 
specified for each mirror. The pre-specified mirror loss 
values represent absorption in the mirror substrates and 
coatings, as well as high-angle scattering due to rough- 
ness finer than the resolution of our grids (for which most 
of the scattered power is lost beyond the finite mirror 
apertures, particularly in the long FP-arms). For each 
individual run in this set, the simulated interferometer 
is completely optimized in the sense of Sec. 1111 CI The 
program input parameters common to all runs are those 
listed in Table ID 

First, however, we note that many tests of the sim- 
ulation program have been performed to ensure its 
validity as a realistic model of a LIGO interferometer. 
These tests include: comparisons with numerical sim- 
ulation results in the literature [H m, to verify ba- 
sic operations like beam propagation and diffractive loss 
from finite mirrors; comparisons against modal analy- 
sis methods |^ for simple interferometer imperfections, 
such as mirror tilts; and comparisons against analytical 
methods 39] for computing the effects of Zernike poly- 
nomial (40j mirror surface deformations. We have imple- 
mented anti-aliasing and energy-conservation procedures 
(cf. Sec. 1111 A)| to make sure that the physics of the inter- 
ferometer is being properly simulated; and we have per- 
formed numerous "common-sense" tests to check that the 
program's output results not only make physical sense, 
but also reasonably reflect the types of interferometer im- 
perfections being modeled. Finally, direct comparisons of 
our simulation results with experimental measurements 
for lar ge-s cale interferometers are now becoming possible 
(e.g., [2l[ 22]), and have thus begun to provide useful 
mutual feedback for both experimental and modeling ef- 
forts. 



A. The Realistically-Deformed Mirror Maps Used 
for these Runs 

The mirror maps used in these simulation runs have 
been derived from two measurements of real optical com- 
ponents: the first one, obtained by LIGO from Hughes- 
Danbury Optical Systems, is a phase map of the reflec- 
tion from the polished surface of the "Calflat" reference 
flat mirror used by the AXAF program (e.g., 'iT]) for 
the calibration of their extremely smooth, high-resolution 
conical mirrors; and the second one is a transmission 
phase map of a trial LIGO mirror substrate obtained 
from Corning. Both of these measurements were of un- 
coated, fused silica substrates. Measurement maps of 
fully-coated mirrors were not available for the runs in 
this study. 

Each of these near-LIGO-quality mirror maps (one sur- 
face reflection map and one substrate transmission map) 
were extrapolated into an array of many maps, so that 
we could place deformed surfaces and substrates upon 
all of the interferometer mirrors simultaneously. One of 



us (Y. H.) used a three-step process for creating new 
mirror deformation maps: first, the surface (or sub- 
strate) source map was Fourier transformed into spatial- 
frequency space; then, the relative phases of its Fourier 
components were randomized; and finally, the applica- 
tion of an inverse Fourier transform created a new mir- 
ror with randomized features, but with the same power 
spectrum of deformations^. This process was performed 
many times, to create 15 new surface maps out of the 
initial Calflat map, and 7 new substrate maps out of 
the initial Corning map. Test runs with various group- 
ings of these mirrors have demonstrated to us that differ- 
ent members of a family of randomized mirrors have (as 
expected) very similar characteristics to one another in 
terms of their effects upon interferometer performance, 
and that swapping one for another has little effect upon 
the output results of the simulation program. 

Further preparation steps for the mirror maps were 
taken to adapt them to the appropriate grid parameters 
and mirror aperture dimensions (cf. Table , and also 
to supply related deformation maps for the beamsplitter, 
given its 45° tilt angle and consequently elliptical aper- 
tures; full details of all such preparations are given in 
[2^ ]. Each of the resulting mirror maps were then tilt- 
removed (using appropriate beam spot sizes) with respect 
to normal-incidence laser beams, as per Sec. 1111 C 61 An 
example of one of these surface deformation maps has 
been shown in Fig. 13 

As a final step, it was necessary to create several fami- 
lies of surface maps with different levels of deformations, 
in order to help the LIGO Project evaluate a range of 
mirror polishing specifications for the procurement of the 
core optics, as well as to make room for the not-fully- 
determined effects of mirror coating deformations. To 
accomplish this, the Calflat-derived surface deformation 
height maps were uniformly multiplied by scale factors 
to generate the new families. The original family of sur- 
faces — which possess RMS deformations of ~.6 nm 
when sampled over their central 8 cm diameters — has 
been labeled "A/ 1800" (with A = Ayag = 1-064 ^m). 
The scaled-up families are labeled A/1200, A/800, and 
A/400, respectively. The mirror substrate maps (possess- 
ing RMS deformations of ~ 1.2 nm over their central 8 cm 
diameters), however, were considered likely to represent 
the best quality of fused silica substrates obtainable for 
the first-generation LIGO interferometers 112] , and were 
not re-scaled; all of the deformed-mirror runs discussed 
below were done with this same family of substrate maps. 
Note that a direct comparison of the substrate RMS value 
with that of the surfaces is not informative, since the sub- 
strates are typically sampled less frequently by the elec- 
tric fields (particularly for the carrier light resonating in 



Although this process does not preserve coherent structures 
— e.g., spikes, etc., which might be systematic effects of mirror 
fabrication — such structures were manually added in by us for 
other tests, enabling us to place limits on their significance. 
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the FP arm cavities), and are therefore less significant 
(for the carrier fields, at least) in their effects. 

Given these families of mirror deformation maps, it 
becomes possible to examine the overall performance of 
a coupled-cavity LIGO interferometer in the presence of 
"realistic" mirror deformations, to comprehensively esti- 
mate its true capabilities. 



B. The Results of the Runs 

This study contains results from five separate simu- 
lation runs: one run with perfectly smooth mirror sub- 
strates and surfaces, and four runs with: (i) deformed 
substrate maps for all of the mirrors, plus, (ii) de- 
formed surface maps for all mirrors from, respectively, 
the A/1800, A/1200, A/800, or A/400 families. For all 
cases other than the "perfect mirrors" run, the transmis- 
sion and (reflective-side) reflection maps for each 2-port 
mirror were constructed from one surface phase map and 
one substrate phase map; the program then derives the 
AR-side reflection map from the other maps via energy 
conservation (i.e., the lossless-mirror Stokes condition, 
cf. Sec. lIII AI . The beamsplitter's two reflection and two 
transmission operators were constructed from one sur- 
face map and two substrate maps, with the remaining 
map derived via its generalized energy conservation for- 
mula. 

The results are summarized in Table ^] Several 
of the quantities described in Sections IIII Gl and IIII Dl 
are included, as well as the true ("absolute") carrier 
and sideband power exiting at the beamsplitter signal 
port. The DG values and pole frequencies of the GW- 
sensitivity curves calculated for each run are also given, 
from which one can construct /isn(/) for each case as 
follows (cf. Eq's. |TT1)-((T^1: 

hsNif) = hsNiO) ■ ^Jl + if/fpoicf ■ (14) 

Some of these quantities have also been re-computed in 
the hypothetical case of an idealized output mode cleaner 
functioning at the signal port, which would act to strip 
away all of the non-TEMgo light (contributing only to 
noise) from the exiting beams, while passing all of the 
TEMoo light (contributing all of the signal and some shot 
noise) through to the output photodetector. A photode- 
tector quantum efficiency of 77 — 0.8 was assumed for 
computing the values of ft.sN(/), F, etc. 

Utilizing this output data, one noteworthy result is 
that several effects of deformed optics have a quadratic 
dependence upon the deformation amplitudes^. This is 



^ We note that the signal-generating sideband power at the beam- 
splitter exit port does not worsen quadratically with RMS defor- 
mation levels; and in fact, TablelTTIshows an increase in exit-port 
sideband power, as the deformation levels get very large. This 
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FIG. 5: The carrier contrast defect, 1 — C, plotted versus RMS 
mirror surface deformations. The dashed line is a quadratic fit 
to the points representing the runs performed with, respec- 
tively: "perfect mirrors", then the A/1800, A/1200, A/800, 
and A/400 mirrors. The horizontal lines are upper limits 
on the contrast defect values allowed for the first-generation 
LIGO and enhanced-LIGO interferometers, respectively. 



as expected, since power scattered out of a Gaussian 
beam by mirror roughness scales like the square of the 
roughness amplitude, even when the deformations are 
spread over a range of spatial frequencies |43|. For ex- 
ample, the effects of imperfect mirrors on the FP arm 
cavity power buildup can be expressed in terms of an 
equivalent "effective mirror loss" that would analytically 
reproduce the same amount of circulating FP-arm power. 
Doing this for each of the runs, we obtain an effective 
loss function that increases quadratically (versus mirror 
surface RMS deformation amplitude) from the baseline 
value of ~50 parts per million of "absorptive" loss that 
is put in by hand for each mirror j24j . Similarly, the 
contrast defect (1 — G), which comes from power coupled 
into non-TEMoo beam modes by (longer spatial wave- 
length) optical imperfections, is also well represented by 
a quadratic fit, as is shown in Fig.jS] Lastly, the optimized 
recycling mirror (power) reflectivity, Ri , should decrease 



counterintuitive behavior is due to great effectiveness of the pa- 
rameter optimizations (cf. Sec. 1111 CI — specifically that of the 
Schnupp length asymmetry, which is forced to a larger value for 
highly deformed mirrors, in order to get the sideband fields out of 
the degraded interferometer as soon as possible. This type of be- 
havior for the sideband fields in the presence of highly deformed 
mirrors, though technically correct, may be unduly optimistic 
when considered in the context of a real LIGO system, in which 
the many parameters are not as easily adjustable — if at all 
adjustable — in the experimental system, as they are in the 
simulation program. 
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TABLE 11: Output results for the series of interferometer simulation runs performed using realistic deformation maps for the 
optical surfaces and substrates. A total (pre-modulation) laser input power of 6 Watts is assumed (except where otherwise 
noted), as well as a photodetector quantum efficiency of »7 = 0.8. 



Quantity 




Values For Specified R,un 




Deformed Surfaces (RMS in wavelengtlis) 


Zero 


Ayag/1800 


Ayag/1200 


Ayag/800 


Ayag/400 


Deformed Substrates (Y/N) 


No 


Yes 


Yes 


Yes 


Yes 


Recycling Mirror Reflectivity" 


98.61% 


98.37% 


98.07% 


97.39% 


93.90% 


Schnupp Lengtli Asymmetry (cm)'^ 


9.0 


12.3 


13.5 


15.9 


24.7 


TEMoo Carrier Power, Recycling Cavity'' 


72.40 


61 54 


51 84 


38 41 


16.37 


TEMoo Carrier Power, Fabry-Perot Arm Avg.' 


4726.7 


4012.0 


3374.3 


2491.7 


1042.4 


TEMoo Carrier Power, Exit Port' 


2.70 X lO""^ 


9.29 X 10"" 


1.94 X 10"^ 


4.56 X 10"^ 


1.89 X lO""* 


Total Carrier Power, Exit Port' 


2.15 X 10"^ 


8.53 x lO"-'* 


1.43 X 10"^ 


2.25 X 10"^ 


3.65 X 10"^ 


Carrier Contrast Defect, 1 — C 


6.02 X 10"^ 


2.82 x 10"* 


5.62 X lO"'' 


1.20 X 10"^ 


4.73 X 10"^ 


TEMoo 1-Sideband Power, Recycling Cavity' 


59.08 


28.32 


25.01 


20.17 


10.66 


TEMoo 1-Sideband Power, Exit Port' 


.9067 


.6745 


.6955 


.7344 


.8196 


Total 1-Sideband Power, Exit Port' 


.9071 


.7590 


.7761 


.8082 


.8795 


GW-Response Pole Frequency, /pole (Hz) 


90.32 


90.38 


90.45 


90.61 


91.45 


Modulation Depth, F ° 


0.279 


0.405 


0.455 


0.501 


0.549 


Absolute Carrier Exit-Port Power (mW) 


12.41 


47.1 


77.2 


118.6 


187.8 


Absolute 2-Sideband Exit-Port Power (mW) 


207.1 


358.9 


458.3 


571.3 


737.7 


DC GW-Sensitivity, hsN{0) ° 


4.79 x lO"^* 


5.76 x 10"^"* 


6.41 X 10"^* 


7.59 X lO"^"" 


1.20 X 10"^^ 


Re-Computed Values Assuming Ideal Output Mode Cleaner: 


Modulation Depth, F " 


0.053 


0.078 


0.093 


0.113 


0.156 


Absolute Carrier Exit-Port Power (mW) 


1.6 X 10"^ 


5.6 x 10"^ 


0.12 


0.27 


1.10 


Absolute 2-Sideband Exit-Port Power (mW) 


7.7 


12.2 


17.9 


28.1 


59.6 


DC GW-Sensitivity, /isn(O) ° 


4.61 X 10"^* 


5.01 X 10"^* 


5.48 X lO"^" 


6.40 X 10"^* 


1.00 X 10"^^ 



"Denotes parameter optimized by program, or during post- 
processing. 

'Denotes quantity normalized to 1 Watt of carrier/sideband exci- 
tation liglit power. 



quadratically from its "perfect mirrors" value, since 1- 
i?i_optim is directly proportional to interferometer losses 
(cf. Sec. lIII Crs)) . and the dominant losses (contrast defect 
loss and high-angle scattering in the FP-arms) are both 
quadratically dependent upon mirror deformation RMS. 
A fit of i?i,optim vs. RMS 24] does indeed bear out this 
expected functional form (though not all the way to a 
surface RMS of zero, since other effects then take over 
— substrate deformations, absorptive losses, etc.). 

The most significant issue to address was whether 
LIGO would be able to perform according to Project 
requirements 44], given mirrors with realistic levels of 
optical deformations. That question is answered in the 
affirmative, as demonstrated in Table ^1 and Fig's. \^ 
|S1 for all cases except that of the worst surfaces simu- 
lated here. First of all, the first-generation LIGO in- 
terferometers are required to have a carrier power gain 
of at least 30 in the PRC, and this target is achieved 
in all runs except for the (ultra-conservative) run per- 
formed with A/400 surface deformations. In addition, 
the first-generation interferometers must have a contrast 
defect of 1 — C < 3 X 10""^, a requirement that is also 
satisfied by all runs other then the A/400 case (and any- 
thing better than ~ A/500 would suffice). Furthermore, 
it has been quoted that the "enhanced" LIGO in- 
terferometers should satisfy the more stringent require- 
ment ofl — C < 1x10^^, which would be achieved 
by three of the five simulation runs here (and anything 
better than ^ A/900 would suffice) — an acceptable 
result, especially considering the likelihood of improved 
mirror quality by the time the enhanced interferometers 



are operational. One caveat, however, is that although 
the contrast defect requirements are met for most of the 
runs, a large amount of total power (several hundred mil- 
liwatts) falls upon the output photodetector in all cases, 
especially for runs with highly deformed mirrors in which 
the optimized sideband modulation depth is large, in or- 
der to help the signal compete against increased shot 
noise. If a signal detection apparatus that could han- 
dle this large amount of power cannot be supplied, then 
an output mode cleaner may be needed; Table Ull shows 
that an output mode cleaner would greatly reduce the 
power that the photodetector must accommodate (while 
also improving the GW-sensitivity by ~15%), as long as 
it operates closely enough to an "idealized" performance, 
as described above. 

The most fundamental requirement to be satisfied is 
that the shot-noise-limited sensitivity curve for an in- 
terferometer with deformed mirrors (cf. Eq's. (|11() . (|14() ') 
should fall within the bounds of the GW- strain- equivalent 
noise envelope requirement, which is the official total- 
noise limit set for the (full 4 km baseline) first-generation 
LIGO interferometers In Fig- HI the computed 

/isn(/) curves for these runs are plotted against the data 
points that define the LIGO requirement envelope, with 
all data given in terms of spectral densities. The seis- 
mic, thermal, and shot-noise-dominated regimes of the 
requirement envelope are apparent in the figure, and 
our /isn(/) curves can be compared to the shot-noise- 
dominated region. We conclude, once again, that all 
runs other than the A/400 case succeed in meeting the 
initial-LIGO requirement. The overall conformity of 
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f (Hz) 

FIG. 6: Comparison of the shot noise curves, /isn(/), com- 
puted for each of the interferometer simulation runs, against 
the official GW-strain-equivalent noise envelope requirement 
for the first-generation LIGO interferometers. 

these results indicates that there is a very clear — and 
very strict, though achievable — quality level for the 
core optics that must be reached in order for the first- 
generation LIGO interferometers to achieve their target 
performances. 



V. IMPACT OF OPTICAL DEFORMATIONS 
UPON LIGO SCIENCE CAPABILITIES 

To place the results of our runs into a scientifically 
relevant perspective, we estimate the effects of opti- 
cal deformations upon LIGO's ability to detect gravita- 
tional waves from anticipated astrophysical events. Fo- 
cusing here upon GW-sources that might be detectable 
by the first-generation LIGO interferometers, and consid- 
ering cases for which an improvement in the shot noise 
limit may significantly increase the number of detection 
events (or enlarge the detectable range of some reason- 
ably well- understood scientific parameter), we arrive at 
two good candidates for study: periodic GW's from non- 
axisymmetric pulsars, and burst GW's from the coales- 
cence of black hole/black hole (BH/BH) binaries. 



A. Non-Axisymmetric Pulsars 



noise curve, /isn(/), that is computed (cf. Eq. (|ll(l . with 
data from Table ITTll for each of the simulation runs. 

These summed, spectral density noise curves must be 
converted into useful signal-to- noise expressions. We fol- 
low the conventions of Thorne , in which each of these 
total- noise curves (/isn(/)) is converted to a dimension- 
less expression (/13/yr), and is compared to the "char- 
acteristic strength" (he) of a GW-source. For periodic 
sources, the condition — h^/yi- means that after co- 
incidence detection in two identical interferometers for 
one-third of a year of integration time, a source with 
strength he can be extracted from the Gaussian noise 
with a confidence level of 90%. 

Averaging over all polarizations and orientations of 
the source on the sky, and treating the GW-frequency 
/ (equal to twice the pulsar's rotation frequency) and 
phase as known, Eq. 52a of Q yields: 

^3/yrl/ w 1.7 • VS • VlO-^Hz X h{f) . (15) 

For a pulsar with gravitational ellipticity e, a rotation- 
axis moment of inertia of 10*^ g • cm^, radiating at fre- 
quency / at a distance r from the earth, and averag- 
ing over orientation angles of the source, we have (from 
Eq. 55 of H): 

/.e|/«7.7xl0--.e.(i^).(^)^ (16) 

In Fig. 13 /13/yr is plotted for all of the simulation 
runs^", and displayed with them are two he curves, each 
representing the locus of possible GW-strengths (plotted 
versus frequency, up to / « 2 kHz) for a pulsar with 
given e and r. We have chosen e = 10^^, which should 
be below the "breaking strain" of neutron star crusts [i^ , 
yet may produce a detectable signal. While this value is 
too large for millisecond pulsars given typical limits on 
their rates of GW-induced spin-down (43, it may not be 
unreasonable for newly- formed pulsars, of which it has 
been estimated 0| that there may be ^^25 such "new" 
pulsars in the galaxy, or one every ~few kpc. 

By setting he = ^3/yr s-t the frequency of peak sensitiv- 
ity (/peak) for a given ft.3/yr curve, one obtains the rough 
estimate that in going from the A/400 case to (or very 
nearly to) the perfect mirrors case, the typical "lookout 
distance" to which such a pulsar is detectable increases 
from ~.57 kpc to 2.0 kpc (while /peak increases to ~200 
Hz). This improvement roughly increases the expected 
number of pulsar detections (in a galactic disk distribu- 
tion) by ^ (2.0/. 57)^ ~ 12, and brings the lookout dis- 
tance for the detection of even a single e = 10^^ pulsar 



To facilitate comparisons against theoretical signal es- 
timates, the output data from the simulation program are 
used to create representative interferometer noise curves. 
To that end, we form the quadrature sum (for each run) 
of three functions: the seismic and thermal noise require- 
ment curves (from Fig.^, added to the particular shot 



Th is f igure differs slightly from the corresponding one (Fig. 2) 
in 14^, because a somewhat older formulation (Eq's. 12 and 13 
of l47l . each with factor of 2 corrections) was used there for the 
thermal noise requirement curve. The numbers quoted in the 
ensuing discussion here differ slightly as a consequence. 
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FIG. 7: Plots of characteristic GW-signal strength he ver- 
sus frequency, for pulsars with specified ellipticity and dis- 
tance from the earth (dashed lines), displayed against the 
dimensionless noise curves /la/yr (for periodic searches) that 
are computed from the output results of the simulation runs 
(solid lines). 

with the first-generation LIGO interferometers to a more 
reasonable value. Alternatively, for a pulsar at a given 
distance and GW-emission frequency, it provides a fac- 
tor of ~ (2.0/. 57) « 3.5 leeway in the smallest detectable 
value of e. 

B. Black Hole/Black Hole Binary Coalescences 

A similar formulation is used for the analysis of burst 
sources. Each quadrature noise sum, h{f), is again com- 
puted; and with appropriate angle averaging and assump- 
tions for optimal filtering, we have (from Eq. 34 of |0|): 

^3/yrl/ « 75 • Vln(//10-« Hz) • V7 X h{f) . (17) 

For bursts, he = h^/yr means that after coincidence de- 
tection in two identical interferometers for one-third of a 
year of observation time, a detection of GW-strength 
has a 90% probability of being a real signal, rather than 
an accidental conspiracy of the Gaussian noise in the de- 
tectors. (Coincidence operation between multiple inter- 
ferometers should theoretically eliminate the false burst 
signals caused by non-Gaussian, uncorrelated noise [26j |. 
which we neglect here.) 

We consider a BH/BH binary system with equal com- 
ponent masses, Mi — M2 — IOMq, located a distance 
r from the earth, and evolving in frequency through the 
inspiral phase until it reaches the onset of the coales- 
cence phase (i.e., merger and ringdown) at /merger ~ 
205 Hz • (20MQ/Mtot) ~ 200 Hz 50]. For these parame- 
ters, Eq. 46b of (3| yields (with a cumulative factor of 2 



adjustment from factor of 2 corrections [HJ to Eq's. 29 
and 44 of (6j|): 

M,.v,.»5.oxio--(MiiP£,,l?^,./..(i„ 

The proper way to interpret this formula, is that if the 
frequency of peak detector sensitivity is / = /peak, then 
the total integrated inspiral signal deposited into the de- 
tector (for comparison with h^/yy) is determined by he 
evaluated specifically at that /peak- 

In Fig. |S1 we have plotted the /13/yr noise curves (for 
burst searches) that are obtained for each of the simula- 
tion runs, along with two he curves (with arrows showing 
time evolution), for coalescence events that just man- 
age to skirt the high-detection-confidence threshold of 
he = ^3/yr a-t /peak during the course of their inspirals. 
The first conclusion that one may draw from this plot 
is that LIGO's sensitivity to these coalescence events (at 
least during inspiral) is most strongly limited by the low 
frequency part of the noise curves, i.e., the seismic and 
thermal noise limits. Nevertheless, there is a measur- 
able benefit from improving the shot noise limit: judging 
from the plot, it can be estimated that in going from the 
A/400 case to the perfect mirrors case, the lookout dis- 
tance is increased from ^--^125 Mpc to ^--^195 Mpc; or equiv- 
alently, it increases the expected number of detectable 
events by the factor ~ (195/125)^ « 3.8. The actual 
rate of BH/BH coalescence events is extremely uncertain 
(even their existence is uncertain), but good middle-of- 
the-road values that one could use as benchmarks are 
the "best estimates" that have been made by Phinney 
[5^. and Narayan et al. [s^, which are, respectively: ~3 
per year out to 200 Mpc (assuming a Hubble constant 
of Hq w 75 km s~^ Mpc~^), and ~1 per year out to 
200 Mpc X (100 km s'^ Mpc~^/Ho). Thus, the perfect 
mirrors run appears to put BH/BH binary coalescence 
events just within the conceivable reach of detection for 
the first-generation LIGO interferometers. Perhaps even 
more significantly, improving the shot noise limit could 
increase LIGO's sensitivity to the onset of the merger 
phase of BH/BH binaries with masses like these; GW- 
emission during the actual merger is still poorly under- 
stood, but it may involve the most powerful radiation of 
detectable energy during the overall coalescence process 
113. 

We end this section by cautioning that the aforemen- 
tioned numbers must be considered very rough estimates, 
given the highly simplified ways in which we have treated 
the noise curves of real interferometers, the sophisticated 
data analysis methods needed to extract signals from the 
noise, and the many inherent uncertainties about the 
GW-sources themselves. In fact, rather than interpret- 
ing the results of this section as indicating how LIGO 
optics would determine initial-LIGO physics, it would be 
more appropriate to interpret these results as a demon- 
stration of how initial-LIGO physics places requirements 
on the optics. The firm scientific conclusion that one can 
draw, however, is that using optical components of the 
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FIG. 8: Plots of characteristic GW-signal strength as a 
function of the detector's peak sensitivity frequency, dur- 
ing the inspiral phase of 10 Mq black hole/black hole bina- 
ries (dashed lines), displayed against the dimensionless noise 
curves /13/yr (for burst searches) that are computed from the 
output results of the simulation runs (solid lines). 

best achievable quality can indeed make a difference in 
whether or not the first-generation LIGO interferometers 
have a fighting chance to detect gravitational waves from 
these promising astrophysical sources. 

VI. DISCUSSION AND CONCLUSIONS 

We have shown our simulation program to be useful 
for gaining physical insight into interferometer behavior, 
and for roughly estimating the effects of optical deforma- 
tions upon LIGO science capabilities. Due to the highly 
detailed nature of the model, it has been an effective tool 
for research and development in the LIGO Project. To 
date, it has been used to address several important de- 
sign issues^^, providing support for technical initiatives 
such as: (i) Aiding LIGO in its transition from Argon-ion 
lasers (A = 514.5 nm) to Nd:YAG lasers (A = 1.064 ^m) 
for the main carrier beam and demonstrating that 
interferometer performance is less sensitive to mirror 
figure deformations (of a specified physical amplitude) 
when larger-wavelength laser beams are used; (ii) As- 
sisting in the selection of the Schnupp Length Asymme- 
try scheme for GW-signal readout over an alternative, 
external modulation (Mach-Zehnder) scheme |0; (iii) 
Helping model the performances (particularly the effects 
of optical deformations upon interferometer control sys- 
tems) of the major LIGO prototype interferometers, in- 



cluding the Fixed Mass Interferometer (FMI) |5^ and 
the Phase Noise Interferometer (PNI) [53 at MIT, and 
the 40-meter interferometer at Caltech (iv) Provid- 
ing assistance in the selection of optical parameters for 
the long-baseline LIGO interferometers, such as beam 
spot sizes, mirror curvatures, and aperture sizes (partic- 
ularly the perspective-foreshortened beamsplitter aper- 
ture) [5^ ; (v) Conducting preliminary studies of the use- 
fulness of an output mode cleaner at the beamsplitter 
signal port (cf. Sec. lIVBll : (vi) Simulating the effects of 
refraction index variations due to thermal lensing (e.g., 
j59|) on interferometer performance, resulting in prelim- 
inary estimates of ^^15% degradation to /isn(/) from 
mirror coating absorption values of 0.6 ppm, or (equiv- 
alently) from substrate bulk absorption values (in fused 
silica) of 5 ppm/cm jBOj . 

Perhaps the most important use of the program has 
been its involvement in the LIGO "Pathfinder Project" 
[6ll |. the initiative to set specifications and tolerances for 
LIGO's core optical components, and to procure them 
through a cooperative effort of several vendors and optics 
metrology groups. Our program has also been used in 
conjunction with other modeling initiatives at LIGO 0, 
to create a broad-based interferometer simulation 
environment involving different algorithmic approaches 
and physical regimes of interest. 

The latest versions of our program continue to be 
used to address important questions raised by the LIGO 
Project, such as estimating the performance of advanced- 
LIGO detectors |^ [53 — i.e., interferometers in- 
corporating Dual Recycling or even Resonant Side- 
band Extraction |^ — in the presence of optical de- 
formations, and participating in initiatives to set core 
optical specifications (and to design the control systems) 
for those advanced detectors [3lll6^ . Many related issues 
will undoubtedly arise in the near future, as advanced in- 
terferometer configurations (and increasingly better op- 
tics) become available, for which this program can be 
used as a primary modeling tool for LIGO and its collab- 
orating gravitational wave groups. 
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